Edit on GitHub

Square Wave Oscillator

Overview

System type: linear hybrid
State dimension: 1
Application domain: Electronics

Model description

Hybrid automaton formulation

using ReachabilityAnalysis, Symbolics, Plots

LazySets.set_ztol(Float64, 1e-15)

function multistable_oscillator(; X0 = Interval(0.0, 0.05),
                                  V₊ = +13.5, V₋ = -13.5,
                                  R = 20.E3, C = 5.5556E-8,
                                  R1 = 20.E3, R2 = 20.E3)

    @variables x
    τ = 1/(R*C)
    α = R2/(R1+R2)

    A = -τ
    b = (τ/α) * V₊
    I₊ = HalfSpace(x <= α*V₊)
    m1 = @system(x' = Ax + b, x ∈ I₊)

    b = (τ/α) * V₋
    I₋ = HalfSpace(x >= α*V₋)
    m2 = @system(x' = Ax + b, x ∈ I₋)

    automaton = LightAutomaton(2)
    add_transition!(automaton, 1, 2, 1)
    g1 = Hyperplane(x == α*V₊)
    r1 = ConstrainedIdentityMap(1, g1)

    add_transition!(automaton, 2, 1, 2)
    g2 = Hyperplane(x == α*V₋)
    r2 = ConstrainedIdentityMap(1, g2)

    modes = [m1, m2]
    resetmaps = [r1, r2]
    H = HybridSystem(automaton, modes, resetmaps)

    # initial condition in mode 1
    X0e = [(1, X0)]
    return IVP(H, X0e)
end
multistable_oscillator (generic function with 1 method)

Results

prob = multistable_oscillator()

sol = solve(prob, T=100e-4, alg=INT(δ=1.E-6), fixpoint_check=false);

plot(sol, vars=(0, 1), xlab="t", ylab="v-")
tspan.(sol)
19-element Array{IntervalArithmetic.Interval{Float64},1}:
 [0, 0.000320001]
    [0.000316999, 0.000888001]
    [0.000883999, 0.00145601]
    [0.00145099, 0.00202401]
    [0.00201799, 0.00259201]
    [0.00258499, 0.00316001]
    [0.00315199, 0.00372801]
    [0.00371899, 0.00429601]
    [0.00428599, 0.00486401]
    [0.00485299, 0.00543201]
    [0.00541999, 0.00600001]
    [0.00598699, 0.00656801]
    [0.00655399, 0.00713601]
    [0.00712099, 0.00770401]
    [0.00768799, 0.00827201]
    [0.00825499, 0.00884001]
    [0.00882199, 0.00940801]
    [0.00938899, 0.00997601]
    [0.00995599, 0.0100211]
location.(sol)
19-element Array{Int64,1}:
 1
 2
 1
 2
 1
 2
 1
 2
 1
 2
 1
 2
 1
 2
 1
 2
 1
 2
 1

Let us analyze in some detail the first transition. If we plot the last 10 reach-sets of the first flowpipe, we observe that only the last 3 actually intersect the guard:

plot(sol[1][end-10:end], vars=(0, 1), xlab="t", ylab="v-")
plot!(x -> 6.75, xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red)

We now cluster those reach-sets into a single hyperrectangle:

Xc = cluster(sol[1], [318, 319, 320], BoxClustering(1))
1-element Array{ReachSet{Float64,Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}},1}:
 ReachSet{Float64,Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}}(Hyperrectangle{Float64,Array{Float64,1},Array{Float64,1}}([6.7258286982437765], [0.024171301756223507]), [0.000316999, 0.000320001])

Plotting Xc matches with the flowpipe after the jump:

plot(sol[1][end-10:end], vars=(0, 1))
plot!(sol[2][1:10], vars=(0, 1))
plot!(x -> 6.75, xlims=(3.1e-4, 3.3e-4), lab="Guard", lw=2.0, color=:red)
plot!(Xc[1], vars=(0, 1), c=:grey)

Finally, we note that the algorithm finds an invariant of the system after the first period. To activate such check pass the fixpoint_check=true flag to the hybrid solve API.

sol = solve(prob, T=100e-4, alg=INT(δ=1.E-6), fixpoint_check=true);
plot(sol, vars=(0, 1), xlab="t", ylab="v-")

When the fixpoint check is activated, the computation terminates as soon as the last reach-set is contained in a previously explored initial state.

tspan(sol)
[0, 0.00145601]

References